clear
version 14.0
set more off
set scheme s2mono


* Specify the year of cash rent data to use in the regression
local reg_year 2013

/*
cd ../maps/
shp2dta using county_my_region2, ///
		database(county_db_my_region2) coordinates(county_coord_my_region2) genid(id) gencentroids(c) replace
*/

cd ../codeAnalysis/
use ../maps/county_db_my_region2, clear
rename stfip stfips
rename cofip cofips
merge 1:1 stfips cofips using "..\dataAnalysis\rent_climate_data", keep(match) nogen
gen ln_cash_rent=ln(cash_rent_non`reg_year')
drop if ln_cash_rent==.
drop if croplandNon_acres==.
spmat contiguity ccounty using ../maps/county_coord_my_region2, id(id) normalize(minmax) replace
spmat idistance dcounty x_c y_c, id(id) normalize(minmax) replace
spmat summarize dcounty

* put spatial weight matrix into a format to be able to calculate Moran's I
spmat export dcounty using "../temp/dcounty.txt", noid replace
insheet using "../temp/dcounty.txt", delim(" ") clear
drop in 1
save "../temp/dcounty.dta", replace
set matsize 1000
spatwmat using "../temp/dcounty.dta", name(dcountyWs) standardize


*capture log close
*log using "..\logs\main results `reg_year'.txt", text replace

* run the do-files to define programs used in this code
run "..\codeSetup\soils transform - program.do"
run "main results - program.do"

local bootstrap_reps 1000

*------------------------------------------------------------------------
* Potential soils
*------------------------------------------------------------------------
local soils  slope om sar gypsum cec7 ec  ksat clay_perc silt_perc ///
			bulkDensity soc0_150 rootznemc ph_less6 ph_greater7dot5 hydgrp* text_*

local soils_nonlin  slope_* om_* sar_* gypsum_* cec7_* ec_*  ksat_* clay_perc_* silt_perc_* ///
			bulkDensity_* soc0_150_* rootznemc_* ph_less6 ph_greater7dot5 hydgrp* text_*

*------------------------------------------------------------------------
* Controls to include in the regression
*------------------------------------------------------------------------
local deficit water_deficit
local surplus water_surplus
local edd dday34C
local soils_controls soc0_150_4 bulkDensity_4 ec_2 ph_less6 ph_greater7dot5 slope_1     
*local soils_controls $SelSoils

local knots_`deficit' 2
local knots_`surplus' 2
local knots_gdd 3
local knots_`edd' 3

local main_year 2013

*------------------------------------------------------------------------
* Regression
*------------------------------------------------------------------------

use "..\dataAnalysis\rent_climate_data", clear
gen ln_cash_rent=ln(cash_rent_non`reg_year')
rename stfips stfip
rename cofips cofip
merge 1:1 stfip cofip using "../maps/county_db_my_region2", keep(match) nogen
rename stfip stfips
rename cofip cofips

capture drop `deficit'_spline*
gen `deficit'_spline = `deficit'
matrix `deficit'_knots=[0]

capture drop water_surplus_spline*
gen `surplus'_spline = `surplus'
matrix `surplus'_knots=[0]

capture drop gdd_spline*
mkspline gdd_spline = gdd, cubic nknots(`knots_gdd') 
matrix gdd_knots=r(knots)

capture drop `edd'_spline*
mkspline `edd'_spline = `edd', cubic nknots(`knots_`edd'') 
matrix `edd'_knots=r(knots)		


summ `deficit'_spline* `surplus'_spline* gdd_spline* `edd'_spline* `soils_controls' stfips_clust croplandNon_acres 

drop if ln_cash_rent==.
drop if croplandNon_acres==.
tab stfips, gen(stfe)
summ `deficit'_spline* `surplus'_spline* gdd_spline* `edd'_spline* `soils_controls'	


reg ln_cash_rent `deficit'_spline* `surplus'_spline* gdd_spline* `edd'_spline* `soils_controls' i.stfips ///
			[aw=croplandNon_acres]
estat hettest
predict ehat, resid
spatgsa ehat, weights(dcountyWs) moran


* Spatial error with heteroskedasticity
* Inverse distance weight matrix
spreg gs2sls ln_cash_rent `deficit'_spline* `surplus'_spline* gdd_spline* `edd'_spline* `soils_controls' stfe1-stfe21 ///
		, id(id) heteroskedastic elmat(dcounty)	
esttab  using "..\tables\spatial_results.tex", ///
			b(%9.3f) se(%9.3f) tex se label replace starlevels(* 0.10 ** 0.05) ///deficit
			r2 sfmt(%9.2f) nogaps mtitles("Coefficients") nonumbers ///
			coeflabels(	water_deficit_spline "Water Deficit" ///
						water_surplus_spline  "Water Surplus"  ///
						gdd_spline1 "GDD Spline 1" gdd_spline2 "GDD Spline 2"  gdd_spline3 "GDD Spline 3"  ///
						dday34C_spline1 "EDD Spline 1" dday34C_spline2 "EDD Spline 2" ///
						soc0_150_4 "SOC" bulkDensity_4 "Bulk Density" ///
						ec_2 "EC" ph_less6 "pH Less than 6" ph_greater7dot5 "pH Greater than 7.5" ///
						slope_1 "Log of Slope") ///
			drop(stfe*)
		
		
* Contiguity weight matrix	
spreg gs2sls ln_cash_rent `deficit'_spline* `surplus'_spline* gdd_spline* `edd'_spline* `soils_controls' stfe* ///
		, id(id) heteroskedastic elmat(ccounty)			

/*
* SARAR	
spreg gs2sls ln_cash_rent `deficit'_spline* `surplus'_spline* gdd_spline* `edd'_spline* `soils_controls' stfe* ///
		, id(id) heteroskedastic dlmat(ccounty) elmat(ccounty)			
			
spreg ml ln_cash_rent `deficit'_spline* `surplus'_spline* gdd_spline* `edd'_spline* `soils_controls' stfe* ///
		, id(id) dlmat(ccounty) elmat(ccounty)						
			
spreg gs2sls ln_cash_rent `deficit'_spline* `surplus'_spline* gdd_spline* `edd'_spline* `soils_controls' stfe* ///
		, id(id) heteroskedastic dlmat(dcounty) elmat(dcounty)	
			
spreg ml ln_cash_rent `deficit'_spline* `surplus'_spline* gdd_spline* `edd'_spline* `soils_controls' stfe* ///
		, id(id) dlmat(dcounty) elmat(dcounty)				
			
*/			
